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Coexisting stochastic and coherence resonance in a mean-field 
dynamo model for Earth's magnetic field reversals 
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PACS 91.25.Mf - Magnetic field reversals: process and timescale 

PACS 91 .25.Cw - Origins and models of magnetic fields; dynamo theories 
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Abstract. - Using a spherical symmetric mean field a 2 -dynamo model for Earth's magnetic 
field reversals, we show the coexistence of the noise-induced phenomena coherence resonance and 
stochastic resonance. Stochastic resonance has been recently invoked to explain the 100 kyr 
periodicity in the distribution of the residence time between reversals. The comparison of the 
resulting residence time distribution with the paleomagnetic one allows for some estimate of the 
effective diffusion time of the Earth's core which may be 100 kyr or slightly below rather than 200 
kyr as it would result from the molecular resistivity. 
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There is ample paleomagnetic evidence that the Earth's 
magnetic field has undergone polarity changes many times 
[1]. Averaged over the last few million years the mean 
rate of reversals is approximately 4-5 per Myr, although 
1 the last reversal occurred approximately 780000 years ago. 
At least two so-called superchrons have been identified as 
periods of some tens of millions of years with no reversal 
, at all. One of the most intriguing features of reversals is 
their pronounced asymmetry with the initial decay of the 
dipole being much slower than the subsequent recreation 
of the dipole with opposite polarity [2]. Despite the gen- 
eral irregularity of the reversal time series, there are (at 
least) two features pointing to some sort of ordering. The 
first one is the clustering property of reversals [3,4], the 
second one is the appearance of a 100 kyr periodicity in the 
distribution of the residence times between reversals [5]. 

Computer simulations of the geodynamo in general, and 
of reversals in particular [6] , have matured much since the 
first fully coupled 3D simulation of a reversal by Glatz- 
maier and Roberts in 1995 [7]. However, the severe prob- 
lem that dynamo simulations have to work in parameter 
regions which are far away from the real values of the 
Earth's core, will remain for a long time. Hence it is cer- 
tainly useful to understand the essential features of rever- 
sals in the framework of simpler models. 

With view on the recent successes of liquid sodium dy- 
namo experiments [8,9] simple reversal models may also 
help in tailoring future dynamo facilities to make them 
prone to reversals. As a matter of fact, reversals were 



observed recently [10] in a special version of the French 
VKS dynamo experiment that was dominated by the use 
of iron propellers with a high magnetic permeability [11]. 
However, the exponential field decay in the initial phase of 
the experimentally observed reversals seems more indica- 
tive for a intermittent extinguishing of the dynamo which 
is certainly not equivalent to the behaviour of the geo- 
dynamo during a reversal. Another interesting dynamo 
related experiment showing reversals has been reported 
by Bourgoin et al. [12]. 

Roughly speaking, there are two classes of simplified dy- 
namo models which try to explain reversals. The first one 
relies on the assumption that the dominant dipole field 
is somehow "rotated" from a given orientation to the re- 
versed one via some intermediate state (or states) [13-15]. 
A typical feature of such models, in particular of the model 
of Hoyng and coworkers [13, 14] , is that they work with an 
effective conductivity of the Earth's core which is strongly 
decreased compared to the molecular value. The necessity 
to introduce such a "turbulent conductivity" is well known 
for the solar dynamo [16], but a dramatic conductivity re- 
duction seems hardly justified for the Earth's core [17]. 
Another drawback of the "rotation model" is that it results 
in reversals with the wrong asymmetry as it was shown 
recently [18]. In spite of these problems, the "rotation 
model" with a turbulent diffusion time of 5 kyr (instead 
of 200 kyr as it would result from the molecular resistiv- 
ity) was successful in recovering the above mentioned 100 
kyr periodicity in the distribution of the residence times 
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between reversals [19]. The underlying physical mecha- 
nism was identified as the stochastic resonance (SR) [20], 
which is an example of noise-induced phenomena in non- 
linear systems driven by weak periodic forcing. For the 
geodynamo, the weak periodic forcing is actually assumed 
to result from the Milankovich cycle of the Earth's orbit 
eccentricity [21,22], although details of the driving mech- 
anism are hard to quantify from first principles. 

The second class of simplified reversal models, which 
could be coined "oscillation models", relies on a spec- 
tral peculiarity of the (in general) non-selfadjoint dynamo 
operator. The basic idea, the specific interplay between 
a non-oscillatory and an oscillatory branch of the domi- 
nant axial dipole mode, was expressed early by Yoshimura 
[23] and later exemplified within a spherically symmetric 
mean-field dynamo model of the a 2 -type [18,24-26]. Many 
features of reversals, in particular the correct time-scale, 
the mentioned asymmetry and clustering property were 
attributed to the magnetic field dynamics in the vicinity 
of a branching point (or exceptional point [27]) of the spec- 
trum of the dynamo operator. Usually, this exceptional 
point, where two real eigenvalues coalesce and continue 
as a complex conjugated pair of eigenvalues, is associated 
with a nearby local maximum of the growth rate situated 
at a slightly lower magnetic Reynolds number. It is the 
negative slope of the growth rate curve between this local 
maximum and the exceptional point that makes station- 
ary dynamos vulnerable to noise. Then, the instantaneous 
eigenvalue is driven towards the exceptional point and be- 
yond into the oscillatory branch where the sign change 
happens. In this picture, reversals appear as noise-induced 
relaxation oscillations. 

In [24] the model was also shown to exhibit coherence 
resonance (CR) [28], which is similar to SR but relies on 
the existence of an intrinsic frequency of the system and 
not on an external periodic forcing. Given the capability 
of this model to account for a number of observed reversal 
features, we will ask in the present paper how the mech- 
anism of SR can be implemented and what we can learn 
from its coexistence with CR. The coexistence of both res- 
onance types has been demonstrated for the Fitz-Hugh- 
Nagumo model with the probability density function of 
the residence times showing a transition from SR to CR 
for an increasing noise strength [29]. 

After delineating our simplified mean-field dynamo 
model (more details and discussions can be found in 
[18,24-26]) we will investigate the coexistence of SR and 
CR. Based on some exploration of the space of governing 
parameters, we will end up with a conjecture on the effec- 
tive diffusive time scale of the core. Actually, this effective 
diffusion time is not well known for two reasons. The first 
one is that the conductivity a of the liquid Fe-Ni-Si alloy 
at the pressure of the Earth's core is hard to determine. 
The most recent estimate [30] is a = 4.71 x 10 5 (£1 m) -1 . 
For the Earth's core of radius R = 3480 km this would 
amount to a magnetic diffusion time := [Iqo~R 2 — 227 
kyr. The second reason is the already mentioned possibil- 



ity that the "molecular" conductivity of the material could 
be reduced due to the turbulent flow. This so-called (3 ef- 
fect [31] is very hard to estimate. Actually, there have been 
claims [32] on the measurement of a few percent (3 effect 
in a turbulent liquid sodium flow with magnetic Reynolds 
numbers Rm up to 8, but neither the Riga nor the Karl- 
sruhe dynamo experiments have shown evidence of any 
significant (3 effect at much larger Rm [8]. However, this 
negative result might not exclude some /3 effect in the geo- 
dynamo which works at still higher Rm and with a higher 
degree of turbulence. Although a very strong reduction 
(by a factor 10 or more) is safely excluded by geomagnetic 
data, a reduction of a by a factor 2 or so is not completely 
forbidden and will be a matter of consideration in this 
paper. 

Our starting point is the well known induction equa- 
tion [31] for evolution of the magnetic field B under the 
influence of a helical turbulence parameter a: 
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After decomposing the magnetic field into a poloidal and 
a toroidal parts according to 



B = — V x (r x VS) - rxVT 



(2) 



the two defining scalars S and T are expanded in spherical 
harmonics of degree I and order m. Under the assumption 
that a is spherically symmetric (which is certainly a grave 
simplification that does not apply to the Earth's outer 
core) we arrive, for each degree I and order m separately, 
at the following pair of partial differential equations: 
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Note that in the non-dimcnsionalizcd equation system (3- 
4) the radius r is measured in units of R, the time r in 
units of the diffusion time Tj, and a in units of (/zoci?) -1 . 
Note further that the order m of the spherical harmonics 
docs not show up in equations (3-4) due to the presup- 
posed spherical symmetry of the problem. The boundary 
conditions are: dsi/dr\ r= i + (I + l).s;(l) = t;(l) = 0. In 
the following we will consider only the dipole field with 
1 = 1. 

Beyond a critical intensity of a, the equation system 
(3-4) will acquire an exponentially growing solution. In 
reality, however, the magnetic field cannot grow indefi- 
nitely, but will attenuate the source of its own generation. 
We model this attenuation by "quenching" the kinematic 
a with the angle averaged magnetic field energy which can 
be expressed in terms of s(r) and t(r). In addition, we as- 
sume that a is also influenced by some noise which might 



p-2 



Stochastic and coherence resonance in a dynamo model for Earth's magnetic field reversals 



co 
c 
CD 
Q 



CO 
X! 
O 




0.4 0.6 
x [Myr] 



Fig. 1: RTD for C = 50, T d = 190 kyr and different values of 
D in the case without periodic driving, compared to the "real" 
data from paleomagnetic measurements which were extracted 
from [5,19]. 



be considered as a shorthand for changing boundary con- 
ditions or the neglected influence of higher magnetic field 
modes. Taking into account quenching and multiplicative 
noise together, we get a time dependent a(r, t) in the form 



a(r, t) 



Oikin{r) 



l + E 



Hi(r) + Ur)r 2 + 6(r)r 3 + £ 4 (r)r 4 



(5) 



with the noise correlation given by 

<&(t)&(t + ti)> = D 2 (l-|n|/T c ) 

x6(l- Inl/TJSij , (6) 

(6 is the Heaviside function). In Eqs. (5,6), a/ci„(r) is 
the kinematic a profile, D is the noise intensity, E is a 
constant measuring the inverse mean field energy, and T c 
is a correlation time of the noise. As already mentioned 
the diffusive time scale Td is approximately 200 kyr if we 
assume the molecular conductivity of the material in the 
Earth's core. However, in the following we will consider 
this value as variable. 

Motivated by the earlier observation that a profiles with 
one sign change along the radius can provide oscillatory 
solutions [33,34], we choose for the kinematic a profile in 
Eq. (5) the Taylor expansion 

ct k in{r) = 1-914 • C ■ {a + a x r + a 2 r 2 + a 3 r 3 + a 4 r 4 ) (7) 

with cto = 1, oei = (X3 = 0, a 2 = —6 and = 5 (the 
factor 1.914 serves just to normalize C in order to make 
it comparable to the case of constant a) . 

For most of the following simulations, the dynamo num- 
ber is set to C = 50 which is highly supercritical com- 
pared to the critical value C cr i t = 6.8 which is specific 
to our particular a(r) profile. The equation system (3-5) 
is time-stepped by means of a standard Adams-Bashforth 
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Fig. 2: Normalized variance (NV) of the RTD for various C 
in dependence on D in the case without periodic driving. The 
minimum at D ~ 9 is a typical signature of CR. 



method with radial grid spacing of 0.02 and time step 
length of 2 xf0~ 5 . The correlation time T c has been set 
to 0.005 which would correspond to f kyr in case that the 
diffusion time is set to 200 kyr. The resulting time se- 
ries show reversal sequences quite similar to those of the 
geodynamo [18,25,26]. As an important characteristic of 
these sequences we will determine the distribution of res- 
idence times between two subsequent reversals which we 
will abbreviate by RTD. 

We start with the case without any periodic forcing of a. 
Figure f shows the RTD of the dynamo model for different 
values of the noise intensity D. The diffusion time was set 
to Td = 190 kyr which is just the double of the 95 kyr 
period which was already mentioned in [22] and which 
was also found in [f 9] as a good fit to the paleomagnetic 
data. For the sake of comparison, we have included the 
curve resulting from paleomagnetic measurement as they 
were published in [5, 19] with its typical SR feature of 
several maxima around multiples of 95 kyr. All of the 
numerically resulting curves exhibit a maximum, followed 
by an exponential decrease of the probability. Starting 
with D = 7, for which the RTD has its maximum at ~ 200 
kyr, a trend of this maximum towards smaller values of r is 
visible. For D = 9 it is located at ^150 kyr and for D = 1 4 
it is at ^f00 kyr. However, for those larger values of D 
we observe a much to steep exponential decay, together 
with an unphysical increase of reversal probability for very 
small values of r. This is a first indication that an assumed 
diffusion time of f 90 kyr might be too large to explain the 
first maximum at ~ 95 kyr observed in the paleomagnetic 
data. 

To illustrate a typical feature of CR we use the normal- 
ized variance NV of the residence time distribution, which 
is defined as follows: 



NV 



s/V^ir) 



(8) 



< T > 

Figure 2 shows the dependence of JVV" on the noise 
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strength D. For all considered values of C, there is a min- 
imum of NV at D ~ 9 which is a clear signature for the 
occurrence of CR [28] . At this value of the noise intensity 
we get an optimal amplification of a residence time cor- 
responding to an inherent timescale of the dynamo. The 
increase of NV for values D > 9 is explained by the fact 
that for large D the dynamo is dominated by the noise 
and very short residence times become more probable as 
it was already visible in the curves for D = 9, 14 in fig. 1. 
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which was also shown to fit better to the paleomagnetic 
data [19]. Actually, the periodic input is implemented 
as an additional periodic variation of the helical turbu- 
lence parameter a with the period Th . Naively, one would 
first apply the periodic forcing to the dynamo number C. 
However, simulations have shown that the dynamo does 
not react significantly to such a homogeneous periodic in- 
put. This insensitivity has to do with the quenching effect 
which reduces the kinematic a to some quenched a profile 
which makes the growth rate of the dynamo close to zero. 
Contrary to such a homogeneous periodic forcing of a, a 
periodic change of the shape of the a(r) profile has a more 
pronounced effect. As a simplest attempt, we have chosen 
to add periodic forcing to the first coefficient cto of the 
Taylor expansion of aki n {r) leaving all other coefficients 
unchanged: 

2tt 

aoO) = 1 + ecos(— • r) (9) 
Jo 

For the period Th we chose always the dominant 95 kyr 
contribution, although the total effect of the Milankovich 
cycles is certainly much more complicated [22]. The pa- 
rameter e measures the relative strength of the periodic 
forcing. 



Fig. 3: RTD for C = 50, T d = 190 kyr, and different values D 
in the case with periodic driving with a period of Tn — 95 kyr 
and forcing strength e = 0.3. 
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Fig. 4: RTD for C = 50, T d = 190 kyr, and different values e 
in the case with periodic driving with a period of Tn = 95 kyr 
and noise intensity D = 7. 

After having discussed the pure CR effect, we will now 
examine the additional occurrence of SR. For this pur- 
pose, a weak periodic forcing of a is introduced. Here, 
"weak" means that there are no reversals in the absence 
of noise (note that a strong periodic force would dominate 
the system completely leading to a peak in the residence 
time distribution located at the timescale of this force). 
Although both additional and multiplicative noise had 
been discussed in [19], we focus here only on the effects of 
multiplicative noise which we consider more physical and 
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Fig. 5: RTD for C = 50, T d = 95 kyr, and different values D 
in the case with periodic driving with a period of Th = 95 kyr 
and forcing strength e = 0.1. 

The results for C = 50, T d = 190 kyr, T n = 95 kyr 
are shown in fig. 3 where we vary D for fixed e = 0.3, 
and in fig. 4 where we vary e for fixed D = 7. We see in 
fig. 3, that for D=7 the residence time distribution shows 
an oscillatory behaviour with a maximum at r = 95 kyr 
and further peaks at integer multiples of Th , which decay 
exponentially. Evidently, the oscillatory behaviour of the 
probability becomes less pronounced for higher values of 
D. It is also visible, that for even higher values of D 
the probability at very small r increases. The effect of 
increasing e (for fixed D = 7) is quite similar to the effect 
of increasing D (for fixed e). However, the increase of the 
probability density for very small r does not show up here. 

Despite some similarity of the curves in figs. 3 and 4 
with the paleomagnetic data, with the given Th = 95 kyr 
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Fig. 6: RTD for C = 50, T d = 95 kyr, and different values e 
in the case with periodic driving with a period of Tn = 95 kyr 
and noise intensity D = 6. 
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Fig. 7: RTD for C = 50, T d = 63.3 kyr and different values D 
in the case with periodic driving with a period of Tq = 95 kyr 
and forcing strength e = O.f . 



and Td — 190 kyr it seems hard to find parameters e and 
D which lead to a real good correspondence. This has 
much to do with the fact that for Td > Tq, which applies 
to our example, the oscillatory character of the probability 
distribution starts only at Td as long as D and e are not 
too large. This behaviour had already been demonstrated 
in [29]. To get the first peak at r = Tq = 95 kyr a 
rather large value of e and/or D is required which, in turn, 
leads to a too steep exponential decrease of the probability 
function. 

For this reason, we will test in the following two addi- 
tional values of the diffusion time, Td — 95 kyr = Tq (figs. 
5 and 6) and T d = 63.3 kyr = 2/3 Tq (figs. 7 and 8). For 
T d = 95 kyr, e = O.f, D = 7 (fig. 5), but also for T d = 95 
kyr, e = 0.15, D = 6 (fig. 6) we observe a quite good 
agreement with the palcomagnctic data, though the ratio 
between the first and the second peak is not exactly the 
same. 

The version with T d = 63.3 kyr, e = 0.1, D = 5.5 (fig. 7) 
looks also quite promising, although the minima between 
the probability peaks seem to be more pronounced than 
in the paleomagnetic data. However, this argument might 
not be so relevant since in the paleomagnetic data some 
stronger "smearing out" of the minima could simply result 
from averaging over a long time period (160 Myr), in which 
some parameters must undergo some changes in order to 
explain the long term variation of the reversal rate (cp. [1], 
p. 184). 

Another point is that in the derivation of the histogram 
of the paleomagnetic data [5], the authors had used a 
moving box technique which automatically leads to some 
smoothing of the curve, ff we do the same moving box 
technique to our data, we arrive at the curves shown in 
fig. 9, which are in reasonable correspondence with the 
paleomagnetic curve. 

To summarize, we have shown that the SR phenomenon 
which seems to lay at the root of the 95 kyr periodicity 
of RTD of the Earth's magnetic field reversal can coexist 



with the CR phenomenon which is typical for our spheri- 
cally symmetric a 2 dynamo model. The interesting ques- 
tion is now if the shape of the residence time distribution 
can serve as a proxy for determining the effective diffu- 
sion time of the Earth's core and hence for the reduction 
of the conductivity due to the f3 effect. Although we had 
shown in [18, 26] that a diffusion time of 200 kyr is well 
compatible with the typical decay and recovery times for 
individual reversals (as long as some high super-criticality 
is assumed) we see now that a diffusion time of f 00 kyr or a 
bit below seems more appropriate when it comes to explain 
the stochastic resonance phenomenon. Finally, this has to 
do with the effect that the first peak of the stochastic res- 
onance at ~1 00 kyr is hardly explainable with a diffusion 
time of ^200 kyr, as long as the strength of the periodic 
forcing and/or the noise are not forbiddingly large. 

Further work will focus on a systematic exploration of 
the space of parameters C, D, Td, and e in order to accom- 
modate simultaneously as much paleomagnetic constraints 
as possible. 
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